#### SECTIONS OF THIS .R FILE: 1) Analyses for the manuscript, 2) Analyses for the appendix.
#### You can search by (e.g.) "Table 1"
#### Questions? szacher@usc.edu  


################ PACKAGES NEEDED
library(haven)
library(tidyverse)
library(srvyr)
library(ggthemes)
library(estimatr)
library(texreg)
library(plm)
library(lmtest)


############## Read in dataset (after converting variables of interest from "wide" to "long" panel data form)
# NOTE: Original dataset available here: https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi%3A10.7910/DVN/TOE8I1 

panel_regs <- read.csv("panel_101214_lagDVs_abortion.csv")
nrow(panel_regs) # 28,212

########################################################### ANALYSES

####### FIGURE 1 (preferences by subgroups)

# PARENT STATUS:
cces_abortionprefs_parent <- panel_regs %>% as_survey(weights = c(weight)) %>% group_by(have_child, year) %>%
  summarize(abortion_prefs_mean = survey_mean(abortion_pref, na.rm = T)) %>% filter(have_child>0) %>% 
  mutate(CATEGORY = ifelse(have_child==1, "Parent",
                           ifelse(have_child==2, "Non-parent", NA)))

# CONSERVATIVE v. NON:
cces_abortionprefs_conservatives <- panel_regs %>% filter(ideology < 6) %>% mutate(cons=ifelse(ideology>3, 1, 0)) %>%
  as_survey(weights = c(weight)) %>% group_by(cons, year) %>%
  summarize(abortion_prefs_mean = survey_mean(abortion_pref, na.rm = T)) %>% 
  mutate(CATEGORY = ifelse(cons==1, "Conservative",
                           ifelse(cons==0, "Liberal & Moderate", NA)))

# PRACTICING CHRISTIANS v. NON: church_att_lag < 4 & religion < 5
cces_abortionprefs_pracchrist <- panel_regs %>% 
  mutate(PracChrist=ifelse(church_att < 4 & religion < 5, 1, 0)) %>%
  as_survey(weights = c(weight)) %>% group_by(PracChrist, year) %>%
  summarize(abortion_prefs_mean = survey_mean(abortion_pref, na.rm = T)) %>% 
  mutate(CATEGORY = ifelse(PracChrist==1, "Practicing Christian",
                           ifelse(PracChrist==0, "Non-practicing/Non-Christian", NA))) %>% filter(PracChrist < 2)

### GENDER
cces_abortionprefs_gender <- panel_regs %>% as_survey(weights = c(weight)) %>% group_by(gender, year) %>%
  summarize(abortion_prefs_mean = survey_mean(abortion_pref, na.rm = T)) %>% 
  mutate(CATEGORY = ifelse(gender==1, "Male",
                           ifelse(gender==2, "Female", NA)))

### SINGLE STATUS:
cces_abortionprefs_single <- panel_regs %>% as_survey(weights = c(weight)) %>% group_by(single, year) %>%
  summarize(abortion_prefs_mean = survey_mean(abortion_pref, na.rm = T)) %>% filter(single >= 0) %>% 
  mutate(CATEGORY = ifelse(single==0, "Married",
                           ifelse(single==1, "Unmarried", NA)))

# plots -- isolated:
typeof(cces_abortionprefs_gender)

plot_abortionprefs_parent <- 
  ggplot() +
  geom_point(data = cces_abortionprefs_parent, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                             color = as.factor(CATEGORY)), size = 3) +
  geom_line(data = cces_abortionprefs_parent, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                            color = as.factor(CATEGORY))) +
  labs(title = "Preferences by Parental Status", x = 'Year',
       y = 'Preference Scale (1: Conservative, 4: Liberal)') +
  scale_y_continuous(breaks = seq(1,4, by = 0.2), limits = c(2.2,3.6)) +
  guides(color=guide_legend(title="Parental status")) +
  scale_colour_grey(end = 0.6) +
  theme_fivethirtyeight() +
  theme(legend.position = "bottom", text = element_text(size = 19), axis.title = element_text())
ggsave(filename = "plot_abortionprefs_parent.jpg", plot = plot_abortionprefs_parent)

plot_abortionprefs_conservatives <- 
  ggplot() +
  geom_point(data = cces_abortionprefs_conservatives, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                                    color = as.factor(CATEGORY)), size = 3) +
  geom_line(data = cces_abortionprefs_conservatives, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                                   color = as.factor(CATEGORY))) +
  labs(title = "Preferences by Ideology", x = 'Year',
       y = 'Preference Scale (1: Conservative, 4: Liberal)') +
  scale_y_continuous(breaks = seq(1,4, by = 0.2), limits = c(2.2,3.6)) +
  guides(color=guide_legend(title="Ideology")) +
  scale_colour_grey(end = 0.7) +
  theme_fivethirtyeight() +
  theme(legend.position = "bottom", text = element_text(size = 19), axis.title = element_text())
ggsave(filename = "plot_abortionprefs_conservatives.jpg", plot = plot_abortionprefs_conservatives)

plot_abortionprefs_pracchrist <- 
  ggplot() +
  geom_point(data = cces_abortionprefs_pracchrist, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                                 color = as.factor(CATEGORY)), size = 3) +
  geom_line(data = cces_abortionprefs_pracchrist, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                                color = as.factor(CATEGORY))) +
  labs(title = "Preferences by Religiosity", x = 'Year',
       y = 'Preference Scale (1: Conservative, 4: Liberal)') +
  scale_y_continuous(breaks = seq(1,4, by = 0.2), limits = c(2.2,3.6)) +
  guides(color=guide_legend(title="Status")) +
  scale_colour_grey(end = 0.7) +
  theme_fivethirtyeight() +
  theme(legend.position = "bottom", text = element_text(size = 19), axis.title = element_text())
ggsave(filename = "plot_abortionprefs_pracchrist.jpg", plot = plot_abortionprefs_pracchrist)

plot_abortionprefs_gender <- 
  ggplot() +
  geom_point(data = cces_abortionprefs_gender, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                             color = as.factor(CATEGORY)), size = 3) +
  geom_line(data = cces_abortionprefs_gender, mapping = aes(x = year, y = abortion_prefs_mean, 
                                                            color = as.factor(CATEGORY))) +
  labs(title = "Preferences by Gender", x = 'Year',
       y = 'Preference Scale (1: Conservative, 4: Liberal)') +
  scale_y_continuous(breaks = seq(1,4, by = 0.2), limits = c(2.2,3.6)) +
  guides(color=guide_legend(title="Gender")) +
  scale_colour_grey(end = 0.6) +
  theme_fivethirtyeight() +
  theme(legend.position = "bottom", text = element_text(size = 20), axis.title = element_text())
ggsave(filename = "plot_abortionprefs_gender.jpg", plot = plot_abortionprefs_gender)

# 

####### TABLE 1 (6 main regressions -- all respondents (2), practicing christians (2), conservatives)

lm_abortionprefs_lagDV <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                      white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                      countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                    data = panel_regs) 

lm_abortionprefs_change <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                       + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                       countycode_lag + relig_import_lag + 
                                       church_att_lag, clusters = caseid, data = panel_regs) 

panel_regs_churchmoreCHRISTIANS <- panel_regs %>% filter(church_att_lag < 4 & religion < 5)

lm_abortionprefs_lagDV_MOREchurchCHRISTIANS <- 
  lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + educ_lag+gender_lag + white_10 + black_10 + latin_10 + 
              faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
            clusters = caseid, data = panel_regs_churchmoreCHRISTIANS)

lm_abortionprefs_change_MOREchurchCHRISTIANS <- 
  lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag+gender_lag + white_10 + black_10 + latin_10 + faminc_lag + 
              pid_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
            clusters = caseid, data = panel_regs_churchmoreCHRISTIANS) 

panel_regs_cons <- panel_regs %>% filter(ideology_lag==4|ideology_lag==5) 

lm_abortionprefs_lagDV_Cons <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 +educ_lag +gender_lag+ 
                                           white_10+ black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                           relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_cons) 

lm_abortionprefs_change_Cons <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag +gender_lag+ white_10+ 
                                            black_10 + latin_10 + 
                                            faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + 
                                            church_att_lag, clusters = caseid, data = panel_regs_cons) 

texreg(list(lm_abortionprefs_lagDV, lm_abortionprefs_change, 
            lm_abortionprefs_lagDV_MOREchurchCHRISTIANS, lm_abortionprefs_change_MOREchurchCHRISTIANS,
            lm_abortionprefs_lagDV_Cons, lm_abortionprefs_change_Cons),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV all", "Chance score all","Lag DV church", "Chance score church",
                              "Lag DV conservative", "Chance score conservative"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

#

###### TABLE 2 (alternative explanations)

table(panel_regs$church_att)
panel_regs_churchatt <- panel_regs %>% filter(church_att < 6)
nrow(panel_regs_churchatt) # 20k

table(panel_regs$ideology)
panel_regs_ideology <- panel_regs  %>% filter(ideology < 6)
nrow(panel_regs_ideology) # 27k

table(panel_regs$pray_freq)
panel_regs_pray_freq <- panel_regs  %>% filter(pray_freq < 7)
nrow(panel_regs_pray_freq) # 23k

table(panel_regs$gun_control)
panel_regs_gun_control <- panel_regs  %>% 
  mutate(gun_control2 = case_when(gun_control==1 ~ 1, # revise because original variable coding is confusing
                                  gun_control==2 ~ 3,
                                  gun_control==3 ~ 2),
         gun_control_lag2 = case_when(gun_control_lag==1 ~ 1,
                                      gun_control_lag==2 ~ 3,
                                      gun_control_lag==3 ~ 2),
         gun_controlCHANGE = gun_control2 - gun_control_lag2)
nrow(panel_regs_gun_control)  # 28k

# urban-rural code change
ALTtest1 <- lm_robust(code_change ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+ black_10 + latin_10 + 
                        faminc_lag + pid_lag + factor(year) + relig_import_lag + church_att_lag, 
                      clusters = caseid, data = panel_regs) 
# county_code_change (i.e., any COUNTY RESIDENCE change)
ALTtest2 <- lm_robust(county_code_change ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+ black_10 + latin_10 + 
                        faminc_lag + pid_lag + 
                        factor(year) + relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs) 
# church_att. (NOTE: DIFFERENT DATA FRAME)
ALTtest3 <- lm_robust(church_attCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+ black_10 + latin_10 + 
                        faminc_lag + pid_lag + 
                        factor(year) + countycode_lag + relig_import_lag, clusters = caseid, 
                      data = panel_regs_churchatt) 
# relig_import.
ALTtest4 <- lm_robust(relig_importCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+ black_10 + latin_10 + 
                        faminc_lag + pid_lag + 
                        factor(year) + countycode_lag + church_att_lag, clusters = caseid, data = panel_regs) 
# PID.
ALTtest5 <- lm_robust(pidCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+ black_10 + latin_10+faminc_lag + 
                        factor(year) + countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                      data = panel_regs) 
# child_healthins.
ALTtest6 <- lm_robust(child_healthinsCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10+black_10 + 
                        latin_10+faminc_lag+ pid_lag + 
                        factor(year) + countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                      data = panel_regs) 
# ideology (NOTE: DIFFERENT DATA FRAME)
ALTtest7 <- lm_robust(ideologyCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 +black_10 + latin_10+ 
                        faminc_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                      data = panel_regs_ideology) 
# gay marriage pref
ALTtest8 <- lm_robust(ban_gaymarCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 +black_10 + latin_10+ 
                        faminc_lag + pid_lag +factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
                      clusters = caseid, data = panel_regs) #
# gun control pref (NOTE: DIFFERENT DATA FRAME)
ALTtest9 <- lm_robust(gun_controlCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 +black_10 + latin_10+ 
                        faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
                      clusters = caseid, data = panel_regs_gun_control) #

# all 9
texreg(list(ALTtest1, ALTtest2, ALTtest3, ALTtest4, ALTtest5, ALTtest7, ALTtest8, ALTtest9, ALTtest6),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Urbal-rural status change", "County change", 
                              "Church att change", "Religious import change", "PID change", "Ideology change",
                              "Gay marriage pref change", "Gun control pref change", "Child health ins change"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

#

###### TABLE 3 (marriage & employment)

panel_regs_unmarmar <- panel_regs %>% group_by(caseid) %>% 
  mutate(cont_unmarried_married = case_when(mean(single_lag)==0|mean(single_lag)==1 ~ 1)) %>% filter(cont_unmarried_married==1)
nrow(panel_regs_unmarmar) #26k

lm_abortionprefs_lagDV_contunmarmar <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                                   white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                   countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                                 data = panel_regs_unmarmar)
lm_abortionprefs_change_contunmarmar <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                                    + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                    countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                                  data = panel_regs_unmarmar)

panel_regs_employ <- panel_regs %>% group_by(caseid) %>% 
  mutate(cont_employ_unemploy = case_when(mean(employ_lag)==0|mean(employ_lag)==1 ~ 1)) %>% filter(cont_employ_unemploy==1)
nrow(panel_regs_employ) #25k
nrow(panel_regs) #28k

lm_abortionprefs_lagDV_contemploy <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                                 white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                 countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                               data = panel_regs_employ)
lm_abortionprefs_change_contemploy <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                                  + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                  countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                                data = panel_regs_employ)

texreg(list(lm_abortionprefs_lagDV_contunmarmar, lm_abortionprefs_change_contunmarmar,
            lm_abortionprefs_lagDV_contemploy, lm_abortionprefs_change_contemploy),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV unmarmar", "Change unmarmar",
                              "Lag DV contemploynonemploy", "Change contemploynonemploy"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

#












########################################################## SUPPLEMENTARY MATERIAL

### Table A1: Non-Effects of Parenthood on Other Outcomes

# (same as above, just with all covariates displayed in the regression table)

# 

### Table A2: Effects of Parenthood among Continuously Unmarried or Continuously Married Respondents.

# (same as above, just with all covariates displayed in the regression table)

# 

### Table A3: Effects for Liberals and Non-Practicing or Non-Christian Respondents.

panel_regs_lib <- panel_regs %>% filter(ideology_lag==1|ideology_lag==2)

lm_abortionprefs_lagDV_Lib <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 +educ_lag + gender_lag+
                                          white_10+ black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                          relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_lib) 

lm_abortionprefs_change_Lib <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag +gender_lag+ white_10+ 
                                           black_10 + latin_10 + 
                                           faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + 
                                           church_att_lag, clusters = caseid, data = panel_regs_lib)

panel_regs_churchless_nonchrist <- panel_regs %>% filter((church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13)) 

lm_abortionprefs_lagDV_LESSchurch <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag+birthyr_10+gender_lag+educ_lag+ 
                                                 white_10+ black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                                 relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_churchless_nonchrist)

lm_abortionprefs_change_LESSchurch <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag+gender_lag+ white_10+ 
                                                  black_10 +  latin_10 + 
                                                  faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + 
                                                  church_att_lag, clusters = caseid, data = panel_regs_churchless_nonchrist) 

texreg(list(lm_abortionprefs_lagDV_Lib, lm_abortionprefs_change_Lib, 
            lm_abortionprefs_lagDV_LESSchurch, lm_abortionprefs_change_LESSchurch),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV liberal", "Change score liberal",
                              "Lag DV more church", "Change score more church"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A4: TWFE regressions

panel_regs_plm_churchless <- panel_regs %>% filter(church_att_lag > 3 & church_att_lag < 7) %>% 
  mutate(new_had_child = ifelse(year>2011, had_child, ifelse(year==2010, 0, NA)))

panel_regs_plm_churchmore <- panel_regs %>% filter(church_att_lag < 4) %>% 
  mutate(new_had_child = ifelse(year>2011, had_child, ifelse(year==2010, 0, NA)))

with(panel_regs_plm_churchless, table(new_had_child, year))
with(panel_regs_plm_churchmore, table(new_had_child, year))
abortion_2wfechurchmore <- plm(abortion_pref ~ new_had_child, 
                               data = panel_regs_plm_churchmore, 
                               index = c("caseid", "year"), model = "within", effect = "twoways") 

abortion_2wfe_regchurchmore <- coeftest(abortion_2wfechurchmore, vcov = vcovHC, type = "HC1", cluster = "group") 

panel_regs_plm_cons <- panel_regs %>% filter(ideology_lag==4|ideology_lag==5) %>% 
  mutate(new_had_child = ifelse(year>2011, had_child, ifelse(year==2010, 0, NA)))

with(panel_regs_plm_cons, table(new_had_child, year))
abortion_2wfe_cons <- plm(abortion_pref ~ new_had_child, 
                          data = panel_regs_plm_cons, 
                          index = c("caseid", "year"), model = "within", effect = "twoways") 

abortion_2wfe_reg_cons <- coeftest(abortion_2wfe_cons, vcov = vcovHC, type = "HC1", cluster = "group") 

texreg(list(abortion_2wfe_regchurchmore, abortion_2wfe_reg_cons),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("2wfe church","2wfe conservative"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A6: Main results among Non-Parents. 

panel_regs_nonparents <- panel_regs %>% filter(havechild_lag==2) # this is non-parents (in t-1)
with(panel_regs_nonparents, table(year, havechild_lag))

lm_abortionprefs_lagDV_nonpar <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                             white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                             countycode_lag + 
                                             relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_nonparents) 
lm_abortionprefs_change_nonpar <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                              + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                              countycode_lag + relig_import_lag + 
                                              church_att_lag, clusters = caseid, data = panel_regs_nonparents) 

texreg(list(lm_abortionprefs_lagDV_nonpar, lm_abortionprefs_change_nonpar),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Results nonparents lag DV","Results nonparents change"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A7a: Main results among Gender Subgroups.

panel_regs_women <- panel_regs %>% filter(gender_lag==2)
panel_regs_men <- panel_regs %>% filter(gender_lag==1)

lm_abortionprefs_lagDV_women <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + 
                                            educ_lag + white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag+ 
                                            relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_women) 

lm_abortionprefs_change_women <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag + 
                                             white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                             relig_import_lag +  church_att_lag, clusters = caseid, data = panel_regs_women) 

lm_abortionprefs_lagDV_men <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + 
                                          educ_lag + white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag+ 
                                          relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_men) 

lm_abortionprefs_change_men <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + educ_lag + 
                                           white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                           relig_import_lag +  church_att_lag, clusters = caseid, data = panel_regs_men) 

texreg(list(lm_abortionprefs_lagDV_women, lm_abortionprefs_change_women,
            lm_abortionprefs_lagDV_men, lm_abortionprefs_change_men),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV women", "Change women",
                              "Lag DV men", "Change men"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A7b: Main results among Gender Subgroups Using an Interaction Term.

lm_abortionprefs_lagDV_gend_int <- lm_robust(abortion_pref ~ had_child + had_child*gender_lag + abortion_pref_lag + birthyr_10 + gender_lag + 
                                           educ_lag + white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                           relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs) # 

lm_abortionprefs_change_gend_int <- lm_robust(abortion_prefCHANGE ~ had_child + had_child*gender_lag + birthyr_10 + gender_lag + educ_lag + 
                                            white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                            relig_import_lag +  church_att_lag, clusters = caseid, data = panel_regs) #

texreg(list(lm_abortionprefs_lagDV_gend_int, lm_abortionprefs_change_gend_int),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV all", "Chance score all"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A8: Main results among Practicing Christians and Conservatives using Interaction Terms.

panel_regs_churchmoreCHRISTIANS01 <- panel_regs %>% 
  mutate(prac_christ_01 = case_when(church_att_lag < 4 & religion < 5 ~ 1,
                                    (church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13) ~ 0),
         cross_pressure = case_when(church_att_lag < 4 & religion < 5 & abortion_pref_lag > 1 ~ 1,
                                    (church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13) | 
                                      (church_att_lag < 4 & religion < 5 & abortion_pref_lag==1) ~ 0)) 

lm_abortionprefs_lagDV_MOREchurchCHRISTIANS_interact <- 
  lm_robust(abortion_pref ~ had_child*prac_christ_01 + abortion_pref_lag + birthyr_10 + educ_lag+gender_lag + white_10 + black_10 +latin_10+ 
              faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
            clusters = caseid, data = panel_regs_churchmoreCHRISTIANS01) 

lm_abortionprefs_change_MOREchurchCHRISTIANS_interact <- 
  lm_robust(abortion_prefCHANGE ~ had_child*prac_christ_01 + birthyr_10 + educ_lag+gender_lag + white_10 + black_10 + latin_10+faminc_lag + 
              pid_lag + factor(year) + countycode_lag + relig_import_lag + church_att_lag, 
            clusters = caseid, data = panel_regs_churchmoreCHRISTIANS01) 

panel_regs_cons2 <- panel_regs %>% mutate(cons = case_when((ideology_lag==4|ideology_lag==5) ~ 1, ideology_lag < 3  ~ 0))

lm_abortionprefs_lagDV_Cons_CROSSint2 <- lm_robust(abortion_pref ~ had_child*cons + abortion_pref_lag + birthyr_10 + educ_lag +
                                                     gender_lag+ 
                                                     white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + countycode_lag + 
                                                     relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_cons2) 

lm_abortionprefs_change_Cons_CROSSint2 <- lm_robust(abortion_prefCHANGE ~ had_child*cons + birthyr_10 + educ_lag + gender_lag + 
                                                      white_10 + black_10 + latin_10 + 
                                                      faminc_lag + pid_lag + factor(year) + countycode_lag + relig_import_lag + 
                                                      church_att_lag, clusters = caseid, data = panel_regs_cons2)

texreg(list(lm_abortionprefs_lagDV_MOREchurchCHRISTIANS_interact, lm_abortionprefs_change_MOREchurchCHRISTIANS_interact,
            lm_abortionprefs_lagDV_Cons_CROSSint2, lm_abortionprefs_change_Cons_CROSSint2),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV Interact PracChris/non", "Change Interact PracChris/non",
                              "Lag DV Interact Conserv/non", "Change Interact Conserv/non"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

#

### Table A9: Abortion Preferences among Practicing Christians and Conservatives.

panel_regs_churchmoreCHRISTIANS <- panel_regs %>% filter(church_att_lag < 4 & religion < 5)
table(panel_regs$abortion_pref_lag)
panel_regs_churchmoreCHRISTIANS01 <- panel_regs %>% 
  mutate(prac_christ_01 = case_when(church_att_lag < 4 & religion < 5 ~ 1,
                                    (church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13) ~ 0),
         cross_pressure = case_when(church_att_lag < 4 & religion < 5 & abortion_pref_lag > 1 ~ 1,
                                    (church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13) | 
                                      (church_att_lag < 4 & religion < 5 & abortion_pref_lag==1) ~ 0)) %>% filter(havechild_lag==2)
nrow(panel_regs) # 28k

# Of 1) practicing Christians, pooled across years, what % answer abortion support Q with 1, 2, 3, 4?
panel_regs_churchmoreCHRISTIANS01_1 <- panel_regs_churchmoreCHRISTIANS01 %>% filter(prac_christ_01==1)
nrow(panel_regs_churchmoreCHRISTIANS01_1) # 9k
table(panel_regs_churchmoreCHRISTIANS01_1$abortion_pref)
2125+ 3800 + 1263 +2185 # 9,373 total ... 
2125/9373 # 22.6%
3800/9373 # 40.5%
1263/9373 # 13.5%
2185/9373 # 23.3%
panel_regs_churchmoreCHRISTIANS01_1_conservatives <- panel_regs_churchmoreCHRISTIANS01_1 %>% filter(ideology_lag==4|ideology_lag==5)
nrow(panel_regs_churchmoreCHRISTIANS01_1_conservatives) # 4520
table(panel_regs_churchmoreCHRISTIANS01_1_conservatives$abortion_pref)
(2256+  482 + 289)/(1459 +2256+  482 + 289) # 67%
panel_regs_churchmoreCHRISTIANS01_1_moderates <- panel_regs_churchmoreCHRISTIANS01_1 %>% filter(ideology_lag==3)
nrow(panel_regs_churchmoreCHRISTIANS01_1_moderates) # 1691
table(panel_regs_churchmoreCHRISTIANS01_1_moderates$abortion_pref)
(505 +346+ 701)/(131+ 505 +346+ 701) # 92%
panel_regs_churchmoreCHRISTIANS01_1_liberals <- panel_regs_churchmoreCHRISTIANS01_1 %>% filter(ideology_lag==1|ideology_lag==2)
nrow(panel_regs_churchmoreCHRISTIANS01_1_liberals) # 1023
table(panel_regs_churchmoreCHRISTIANS01_1_liberals$abortion_pref)
(125+ 140+ 719)/(28 +125+ 140+ 719) # 97%

panel_regs_conservatives <- panel_regs %>% filter(ideology_lag==4|ideology_lag==5)
table(panel_regs_conservatives$abortion_pref)
2422 +5337 +1970+ 2044 # 11,773 total ... 
2422/11773 # 20.6%
5337/11773 # 45.3%
1970/11773 # 16.7%
2044/11773 # 17.4%

# 

### Table A10: Main Results including State Fixed Effects.

lm_abortionprefs_lagDV_stateFE <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                      white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                      countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                    fixed_effects = ~state_lag, data = panel_regs) 

lm_abortionprefs_change_stateFE <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                       + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                       countycode_lag + relig_import_lag + 
                                       church_att_lag, fixed_effects = ~state_lag, clusters = caseid, data = panel_regs) 

texreg(list(lm_abortionprefs_lagDV_stateFE, lm_abortionprefs_change_stateFE),
       stars = c(0.01, 0.05),
       custom.model.names = c("Lag DV", "Change"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A11: Descriptive statistics

# DV (1-4)
table(panel_regs$abortion_pref)
2787 + 7335 + 3934 +13962
mean(panel_regs$abortion_pref, na.rm = TRUE) 
sd(panel_regs$abortion_pref, na.rm = TRUE) 
# ideology (1-5)
table(panel_regs$ideology)
2925 +5109+ 7789 +7101 +4588
panel_101214_lagDvs_abortion_ideology <- panel_regs %>% filter(ideology < 6)
mean(panel_101214_lagDvs_abortion_ideology$ideology, na.rm = TRUE)
sd(panel_101214_lagDvs_abortion_ideology$ideology, na.rm = TRUE)
# practicing Christians (0/1)
panel_101214_lagDvs_abortion_pracChrist <- panel_regs %>% 
  mutate(prac_christ_01 = case_when(church_att_lag < 4 & religion < 5 ~ 1,
                                    (church_att_lag > 3 & church_att_lag < 7) | (religion > 4 & religion < 13) ~ 0))
table(panel_101214_lagDvs_abortion_pracChrist$prac_christ_01)
18644 + 9448
mean(panel_101214_lagDvs_abortion_pracChrist$prac_christ_01, na.rm = TRUE) 
# single? (0/1)
table(panel_regs$single)
19533 + 8678 
mean(panel_regs$single, na.rm = TRUE) 
# employ? (0/1)
table(panel_regs$employ_lag)
13749+ 14463
mean(panel_regs$employ_lag, na.rm = TRUE) 
# gender (0/1)
table(panel_regs$gender)
15732+ 12480 
mean(panel_regs$gender, na.rm = TRUE) 
# white (0/1) / black / latin / PID
table(panel_regs$pid_lag)
13015 + 2928 +12203
12203/28145
mean(panel_regs$latin_10, na.rm = TRUE) 
# income 
table(panel_regs$faminc_lag)
panel_101214_lagDvs_abortion_faminc_lag <- panel_regs %>% filter(faminc_lag < 8)
nrow(panel_101214_lagDvs_abortion_faminc_lag)
mean(panel_101214_lagDvs_abortion_faminc_lag$faminc_lag, na.rm = TRUE) 
sd(panel_101214_lagDvs_abortion_faminc_lag$faminc_lag, na.rm = TRUE) 
# birth year
table(panel_regs$birthyr_10)
panel_101214_lagDvs_abortion_birthyr_10 <- panel_regs %>% filter(birthyr_10 < 1993)
nrow(panel_101214_lagDvs_abortion_birthyr_10)
mean(panel_101214_lagDvs_abortion_birthyr_10$birthyr_10, na.rm = TRUE) 
sd(panel_101214_lagDvs_abortion_birthyr_10$birthyr_10, na.rm = TRUE) 
# education
table(panel_regs$educ_lag)
panel_101214_lagDvs_abortion_educ_lag <- panel_regs %>% filter(educ_lag < 7)
nrow(panel_101214_lagDvs_abortion_educ_lag)
mean(panel_101214_lagDvs_abortion_educ_lag$educ_lag, na.rm = TRUE) 
sd(panel_101214_lagDvs_abortion_educ_lag$educ_lag, na.rm = TRUE) 
# urbal-rural
table(panel_regs$church_att_lag)
panel_101214_lagDvs_abortion_countycode_lag <- panel_regs %>% filter(church_att_lag < 8)
nrow(panel_101214_lagDvs_abortion_countycode_lag)
mean(panel_101214_lagDvs_abortion_countycode_lag$church_att_lag, na.rm = TRUE) 
sd(panel_101214_lagDvs_abortion_countycode_lag$church_att_lag, na.rm = TRUE) 

# 

### Table A12: Main Results among Respondents whose Family Income Did Not Meaningfully Change.

with(panel_regs, table(year, faminc_lag))
sd(panel_regs$faminc_lag, na.rm = TRUE) # 1.8
panel_regs_inc <- panel_regs %>% group_by(caseid) %>% 
  mutate(cont_faminc_lag = case_when(sd(faminc_lag)<.45 ~ 1)) %>% filter(cont_faminc_lag==1)
nrow(panel_regs_inc) # 15.3k
lm_abortionprefs_lagDV_continc <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                              white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                              countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_inc)
lm_abortionprefs_change_continc <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                               + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                               countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_inc)

texreg(list(lm_abortionprefs_lagDV_continc, lm_abortionprefs_change_continc),
       stars = c(0.01, 0.05),
       custom.model.names = c("Lag DV continc", "Change continc"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A13: Non-effects on Having a Second+ Child on Abortion Preferences

panel_regs_par <- panel_regs %>% filter(havechild_lag==1)

lm_abortionprefs_lagDV_2ndchild <- lm_robust(abortion_pref ~ had_2nd_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                               white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                               countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, 
                                             data = panel_regs_par) 
lm_abortionprefs_change_2ndchild <- lm_robust(abortion_prefCHANGE ~ had_2nd_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                                + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                countycode_lag + relig_import_lag + 
                                                church_att_lag, clusters = caseid, data = panel_regs_par) 
texreg(list(lm_abortionprefs_lagDV_2ndchild, lm_abortionprefs_change_2ndchild),
       stars = c(0.01, 0.05),
       custom.model.names = c("Lag DV 2nd child all", "Change 2nd child all"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A14: Main Results including State Citizen Ideology as Control.
# Data from https://rcfording.wordpress.com/state-ideology-data/.

state_ideology_2012 <- read.csv("state_ideology_2012.csv")
state_ideology_2012_B <- state_ideology_2012 %>% mutate(state_lag = FIPS) %>% select(state_lag, year, IDEO_1_citizen)
panel_regs_WITH_state_ideology_2012_B <- panel_regs %>% left_join(state_ideology_2012_B, by = c("state_lag", "year"))
table(panel_regs_WITH_state_ideology_2012_B$IDEO_1_citizen)

lm_abortionprefs_lagDV_stateideo <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                                white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                countycode_lag + relig_import_lag + church_att_lag + IDEO_1_citizen, clusters = caseid, 
                                              fixed_effects = ~state_lag, data = panel_regs_WITH_state_ideology_2012_B) 
lm_abortionprefs_change_stateideo <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                                 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                 countycode_lag + relig_import_lag + church_att_lag + IDEO_1_citizen, 
                                               fixed_effects = ~state_lag, clusters = caseid, data = panel_regs_WITH_state_ideology_2012_B) 
texreg(list(lm_abortionprefs_lagDV_stateideo, lm_abortionprefs_change_stateideo),
       stars = c(0.01, 0.05),
       custom.model.names = c("Lag DV state ideo", "Change state ideo"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs are clustered by individual and robust.", 
       include.ci = FALSE) 

# 

### Table A15: Main Results among Married versus Unmarried Respondents

panel_regs_MAR <- panel_regs %>% filter(single_lag==0)
panel_regs_UNMAR <- panel_regs %>% filter(single_lag==1)

lm_abortionprefs_lagDV_contMAR <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                              white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                              countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_MAR)
lm_abortionprefs_change_contMAR <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                               + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                               countycode_lag + relig_import_lag + church_att_lag, clusters = caseid, data = panel_regs_MAR)
lm_abortionprefs_lagDV_contUNMAR <- lm_robust(abortion_pref ~ had_child + abortion_pref_lag + birthyr_10 + gender_lag + educ_lag + 
                                                white_10 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                countycode_lag + relig_import_lag + church_att_lag, clusters = caseid,data = panel_regs_UNMAR)
lm_abortionprefs_change_contUNMAR <- lm_robust(abortion_prefCHANGE ~ had_child + birthyr_10 + gender_lag + educ_lag + white_10 + 
                                                 + black_10 + latin_10 + faminc_lag + pid_lag + factor(year) + 
                                                 countycode_lag + relig_import_lag + church_att_lag, clusters = caseid,data = panel_regs_UNMAR)

texreg(list(lm_abortionprefs_lagDV_contMAR, lm_abortionprefs_change_contMAR,
            lm_abortionprefs_lagDV_contUNMAR, lm_abortionprefs_change_contUNMAR),
       stars = c(0.01, 0.05, 0.1),
       custom.model.names = c("Lag DV mar", "Change mar",
                              "Lag DV unmar", "Change unmar"),
       digits = 2, 
       caption = "Having children & Abortion policy prefs (CCES). SEs clustered by individual and robust.", 
       include.ci = FALSE) 

#











